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ABSTRACT 

We study the problem of time-dependent photoionization of low density gaseous nebulae subjected to sudden 
changes in the intensity of ionizing radiation. To this end, we write a computer code that solves the full time- 
dependent energy balance, ionization balance, and radiation transfer equations in a self-consistent fashion for a 
simplified pure hydrogen case. It is shown that changes in the ionizing radiation yield ionization/thermal fronts that 
propagate through the cloud, but the propagation times and response times to such fronts vary widely and nonlinearly 
from the illuminated face of the cloud to the ionization front (IF). IF /thermal fronts are often supersonic, and in slabs 
initially in pressure equilibrium such fronts yield large pressure imbalances that are likely to produce important 
dynamical effects in the cloud. Further, we studied the case of periodic variations in the ionizing flux. It is found that 
the physical conditions of the plasma have complex behaviors that differ from any steady-state solution. Moreover, 
even the time average of ionization and temperature is different from any steady-state case. This time average is 
characterized by overionization and a broader IF with respect to the steady-state solution for a mean value of the 
radiation flux. Around the time average of physical conditions there is a large dispersion in instantaneous conditions, 
particularly across the IF, which increases with the period of radiation flux variations. Moreover, the variations in 
physical conditions are asynchronous along the slab due to the combination of nonlinear propagation times for 
thermal fronts/IFs and equilibration times. 
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1. INTRODUCTION 

The general problem of photoionization modeling has broad 
importance in astrophysics. This topic comprises any situation 
in which energy in the form of electromagnetic radiation is 
provided to a gaseous object. The radiation is then re-processed 
by the gas, which becomes ionized and heated, and the excess 
energy is re-emitted into longer wavelength spectral lines and 
diffuse continuum. 

Traditionally, modeling of astronomical photoionized plas- 
mas is done from the condition of steady-state statistical 
equilibrium, which means that gas ionization is balanced by 
recombination, atomic excitations are balanced by spontaneous 
and induced de-excitations, and electron heating is balanced 
by cooling. These conditions result in coupled ionization/ 
excitation balance equations (one for each atom and ion in 
the plasma) and a general thermal balance equation. In addi- 
tion, the models must determine the local radiation field, in- 
cluding direct and diffuse components, which is also coupled 
to the conditions above through the radiative transfer equation 
(Osterbrock & Ferland 2006). There has been much progress in 
steady-state photoionization modeling in the last few decades 
through increasingly detailed treatment of the microphysics, 
improvements in the quality and completeness of atomic and 
molecular data, and growth of computational power. At present 
there are several sophisticated photoionization modeling codes 
in use, e.g„ XSTAR (Kallman & Bautista 2001), CLOUDY 
(Ferland et al. 1998), TLUSTY (Hubeny & Lanz 1995), and 
MOCASSIN (Ercolano et al. 2003). 

3 Current address: Harvard-Smithsonian Center for Astrophysics, 60 Garden 
Street, Cambridge, MA 02138, USA. 


The steady-state assumption is appropriate whenever the 
equilibration timescales for excitation, ionization, and thermal 
balance are much shorter than variability timescales in either 
the ionizing radiation continuum or the geometrical structure 
of the plasma. However, if the ionizing radiation changes at a 
rate slower than the equilibrium timescales, or if other condi- 
tions change on timescales shorter than those of microscopic 
processes, then it is necessary to take into account the full 
temporal dependence of the state equations. There are many 
astrophysical systems in which time-dependent photoioniza- 
tion (TDP) modeling has been discussed. Some examples in- 
clude the interstellar medium (Lyu & Bruhweiler 1996; Joulain 
et al. 1998), Hu regions (Rodriguez-Gaspar & Tenorio-Tagle 
1998; Richling & Yorke 2000), planetary nebulae (Harrington & 
Marionni 1976; Harrington 1977; Schmidt- Voigt & Koeppen 
1987; Frank & Mellema 1994; Marten & Szczerba 1997), novae 
and supernovae (Hauschildt et al. 1992; Beck et al. 1995; Kozma 
& Fransson 1998; Dessart & Hillier 2008), reionization of the 
intergalactic medium (Ikeuchi & Ostriker 1986; Shapiro & Kang 
1987; Shapiro et al. 1994; Ferrara & Giallongo 1996; Giroux & 
Shapiro 1996; Ricotti et al. 2001), ionization of the solar chro- 
mosphere (Carlsson & Stein 2002), gamma ray bursts (Perna 
& Loeb 1998; Bottcher et al. 1999), accretion disks (Woods 
et al. 1996), active galactic nuclei (Nicastro et al. 1997, 1999; 
Krongold et al. 2007), evolution of the early universe (Seager 
et al. 2011), and quasar FeLoBALs (Bautista & Dunn 2010). 
However, there is as yet no general tool to model non- 
equilibrium photoionized plasmas. 

In this paper we lay out the basic approach to solve the 
TDP problem and present an overview of the behavior of non- 
equilibrium, pure hydrogen photoionized plasmas. We illustrate 
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the behavior in various cases of general interest in astrophysics. 
This is a first step toward the development of a general purpose 
TDP modeling code. 

While the treatment of pure hydrogen plasmas does not in- 
clude all the complexity of chemically enriched nebulae it is 
interesting to study this case in detail. First, the treatment of 
time-dependent effects, while expected to be present in vari- 
ous scenarios, implies introducing a number of new parameters 
and complexities for nebular modeling. Thus, it is important 
to introduce time-dependent effects progressively in order to 
understand the physics in detail and being able to disentan- 
gle time-dependent effects from already known variables like 
optical depth, adopted spectral energy distributions, chemical 
effects, etc. Thus, an extensive study of optically thin, pure hy- 
drogen nebulae, however qualitative, is a natural and necessary 
first step in toward time-dependent modeling. Another motiva- 
tion for this study is the ongoing Z-pinch experiments, like those 
at the University of Nevada (e.g., Mancini 2011), which seek to 
test the accuracy of photoionization modeling codes on single 
composition plasmas, for example pure hydrogen or pure neon, 
or particular mixtures of two gases. However, these experiments 
are intrinsically time-dependent. 

2. FUNDAMENTAL EQUATIONS 
2.1. Ionization Balance 


In this simplified model there is one free electron per bare 
proton, i.e., n 2 = n e . Furthermore, given that the hydrogen 
density, n = n\ + n 2 , is conserved Equation (1) can be written 
in terms of ;?i as 

= n 2 ( x , t) [a r (T) + a c (T)] - m(x, t){y(x , t ) 
at 

+ n[2a r (T) + a c (T)]} + n 2 a r (T) (5) 

2.2. Energy Equation 

The temperature of the gas is found by solving the energy 
equation. The net heat of the system is given by 

dQ _ ^(heat) _ p(cool) /g\ 

dt 

where Q is the particle kinetic energy and the terms on the 
right-hand side of the equation are the heating and cooling rates. 
Here, we consider heating by photoionization and cooling by 
recombination and collisional ionization. 

By assuming rapid energy equipartition among atoms, pro- 
tons, and electrons, one can write the particle kinetic energy as 
Q = ( 3/2)n t kT , where n t = n + n e = 2n — n\ is the total 
number density, k is the Boltzmann constant, and T is the gas 
temperature. Then, if the number density, n, is constant one finds 


As a first approach to an otherwise cumbersome problem, we 
will start by considering a gas composed entirely of hydrogen. 
Additionally, we approximate the system to only two energy 
levels, i.e., one to represent the ground state and another to 
represent the continuum. This means that no bound excited states 
are included, and only ionization and recombination processes 
are considered. Under these assumptions, the population n\ of 
the ground state can be described as 


dn i(x, t) 
dt 


—n i(jc, t) [y(x, t ) + n e a c (T)\ + n 2 (x, t)n e a r (T), 


( 1 ) 

where n, is the population of level i and n e is the electron density. 
y is the photoionization rate, which is given by 


y(x, t ) 


-f 


’ ds 

( r e J E (x , t ) — , 
e 


( 2 ) 


where er E is the photoionization cross section and J s (x, t ) is the 
mean intensity of the radiation field. a c and a r are the collisional 
ionization and recombination rate coefficients, respectively. For 
the collisional ionization rate, we will adopt the expression given 
in Cen (1992): 


a c (r) = 5.83 x 10“ 11 7’ 1/2 (1 + 7’ 5 1/2 ) V 57809 - i/7 ' (3) 


dT 

dt 


2 

3(2 n — n\)k 


^(heat) _ p(cool) + 

2 


dn\ 

dt 


(7) 


The last term on the right-hand side of this equation corresponds 
to changes in the kinetic energy associated with temporal 
changes in the ionization of the plasma. This term explicitly 
couples the ionization and thermal balance equations, but it is 
zero under steady-state conditions. The photoionization heating 
is 

C ^ d s 

A (ph0) = / o s J (x, t) e n\(x, f)(e — e th ) — (8) 

Jo £ 

and can be written as 

A (pho) = «i(x, t)y(x, t)(e) (9) 


where 


(£> 


/“ JJx, t)a e {e - e lb )de/e 
/“ J E (x, t)a s ds/e 


( 10 ) 


is the mean kinetic energy of free electrons weighted by the pho- 
toionization cross section, and e t h = 13.6 eV is the threshold 
energy for hydrogen. The recombination and collisional ioniza- 
tion cooling rates are given by 


r (rec) = n e n 2 (x, t)a,.(T)gkT (11) 


and for the recombination rate we will use the fitting formula anc j 
given by Badnell (2006), 


T (col) = n e ni(x,t)a c (T)s th . (12) 


a r (T) = 8.32 x 1 0~ n |y 772.97(1 + yr/2.97) 1 - 075 

x (1 -y/T/1 x io 5 ) 1+0 - 75 r\ (4) 


In Equation (1 1) g is a constant factor, typically about 0.6, that 
depends on the spectral energy distribution of the radiation field. 
Then, the thermal balance equation can be written as 


where T is the temperature in Kelvin, T$ is in units of 10 5 K, and 
a r and a c are both in cm 3 s _1 . Note that we only consider the so 
called ‘Case A’ recombination case, which is consistent with the 
assumption of optically thing nebulae. Other scenarios, such as 
Case B recombination of hydrogen, will be treated elsewhere. 


dT(x, t) 
dt 


3(2 n — n\)k 


n\(x, t)y(x, t)(e) — kTn 2 e a r {T) 

03) 


3 dn 1 

-n\{x, t)n e a c (T)e th + -kT — — 
2 dt 
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Figure 1. Time-dependent simulation for a slab with constant density of n = 10 4 cm ~ 3 and initial flux of F x = 7.95 erg cm ~ 2 s -1 . At t = 0 s the flux is increased 
by a factor of three. The upper and lower panels show the neutral hydrogen density and the gas temperature along the position within the slab, respectively. In 
both cases, each curve corresponds to the profile at a different moment in time. The initial condition is plotted in red, and the final state of the system is plotted in 
green. 

(A color version of this figure is available in the online journal.) 


This equation has no analytic solution, even in the steady state 
case dT /dt = 0 due to the nonlinear dependence of a r and a, 
on T. 

In general. Equations (5) and (13) need to be solved simul- 
taneously. Moreover, they both depend on the radiation field, 
which needs to be known at each position and instant in time. 
Thus, one also needs to solve a coupled equation for radiation 
transfer. 


2.3. Ionization Parameter and Radiative Transfer 

For the sake of clarity, it is assumed that the spectral energy 
distribution of the source remains constant. Then, as shown by 
Tarter et al. (1969), the state of the gas is determined by a single 
parameter known as the ionization parameter 

S = * 4ttF x (s), (14) 

nR z 

where (e) is the mean photon energy and R is the distance from 
the source, L is the luminosity (in energy units) of the ionizing 
source, and F x is the flux of ionizing radiation. In practice L 
is integrated from 1 Ry, the ionization threshold for hydrogen, 
to 1000 Ry, beyond which the radiation is expected to be very 


small. L and F x are related through 

1 L v 

F x = / —dv. 

47 rR~ J 1Ry hv 

This definition for the ionization parameter is related to various 
other customary ionization parameter definitions, i.e., IJ H = 
F x /n (Davidson & Netzer 1979); Z = F v {vj/)/{2hcn), where 
F v (Vf ) is incident (energy) flux at 1 Ry; and 5 = L/(4R2cnkT) 
(Krolik et al. 1981). 

The radiation transfer equation describes the interaction of 
the radiation from the source and the material in the gas. In 
plane-parallel geometry, the time-dependent radiative transfer 
equation can be written as 


1 3 I e (x, /r, t) 
c dt 


dI E (x, fi, t) 


+ F 


dx 


r) e (x, t) - Xb(x, t)I e {x, fi, t ), 


(15) 

where I e (x, /r, t) is the intensity of the radiation field, /i is the 
cosine of the angle with respect to the normal, and r) e {x, t ) 
and Xf (x , t ) are the total emissivity and opacity, respectively. 
The solution of this equation is computationally challenging, as 
discussed extensively in the literature. For the present qualitative 
TDP study we simplify this equation by adopting a one-stream 
approximation, in which only the direction along the normal is 
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Figure 2. Ionization (solid line) and recombination (dotted line) rates vs. depth inside the slab with log £ = 0 after a sudden increase of the ionizing flux by a factor 
of three. The rates are plotted at t = 0 (initial steady-state conditions), t = 3.4 x 10 8 s (when the slab has reached equilibrium again), and two instants in between. 
(A color version of this figure is available in the online journal.) 


considered (i.e., /x = 1), and then J, = (1/2) I s dfi ~ l t . 
Furthermore, by neglecting any local emissivity within the gas 
as well as photon scatterings, the radiative transfer equation can 
now be written as 


1 3 J E (x, t ) 
c 3 1 


+ 


3 J E (x, t ) 
dx 


-n i(*, t)a E J E {x , t). 


(16) 


2.4. Characteristic Times 


Note that these definitions of ionization and recombination 
times are different from more conventional definitions in that 
we include the factors n/n\ and n/ni. For example, a typical 
definition of recombination time is f rec = 1 /(n e a r ), which is 
appropriate for the steady-state condition in the fully ionized 
region where n/n 2 ~ 1. Our present definitions are generally 
correct for nebulae where the ionization of the plasma may 
change with time. 

The ionization equilibration time, r; on , can be defined as 


The response of a plasma to variations in an ionizing radiation 
source is governed by three timescales: the ionization equilibra- 
tion timescale, the temperature equilibration timescale, and the 
propagation timescale. 

In terms of the ionization of the plasma we have the photoion- 
ization time 

n 

(17) 


n 

n\Y ’ 


the recombination time 


n 


tree — 

n 2 n e <x r 

and the collisional ionization time 


(18) 


1 (hi - rtf) 

Tl Tj on 


1 1 

Aon free 


(20) 


where nf is the equilibrium neutral hydrogen density after the 
change in radiation field and f Kln is defined as the ionization 
time, tj on — tpi ? C oi /17p: tcoi )■ Thus, 


«r “ ”i 


Ann tre 


( 21 ) 


In terms of the temperature behavior, it is useful to define the 
temperature equilibration time, tj, as 


n 

tcol — • 

n\n e a c 


(19) 


3 k(T - T e ) 

2,Xf 


d 

dt 



( 22 ) 
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Figure 3. Heating (solid line) and cooling (dotted line) rates vs. depth inside the slab with log £ = 0 after a sudden increase of the ionizing flux by a factor of three. 
The rates are plotted at t = 0 (initial steady-state conditions), t = 3.4 x 10 8 s (when the slab has reached equilibrium again), and two instants in between. 

(A color version of this figure is available in the online journal.) 


where T E is the equilibrium temperature after the change in 
radiation field. For constant n and n \ , r T is the ratio of the 
excess of energy density to the net cooling rate A — T. 

The ionization and temperature timescales are intrinsically 
related through various rate coefficients involved. Nonetheless, 
the former is generally much longer than the latter, as illustrated 
in the next section. 

The equilibration timescales defined above refer to changes 
in local conditions under variations in the local radiative field. 
However, such changes are not simultaneous across the cloud. 
Instead, variations in the local radiation field at any depth inside 
the cloud are delayed with respect to the illuminated face of the 
cloud by the radiation propagation time 


‘•pro 


-L 


x n dr) 
F(r) 


dr 


x(n i) 

F x 


Nh 

F x ’ 


(23) 


where Nh is the neutral hydrogen column density; see also 
Schwarz et al. (1972). The propagation time is the characteristic 
time it takes for the ionization front (IF) to move under the 
assumption that there is one ionization event per incident photon. 
The above equation shows that variations in the radiation field 
propagate quickly and at nearly constant rate through the ionized 
region, but the propagation time increases steeply across the IF, 
where n\ increases. Thus, large departures from equilibrium 
conditions should be expected across the IF under variations 
of the radiation field. Across the IF too the equilibration times 


reach maximum values. Thus, the IF expected to exhibit the 
largest departures from equilibrium conditions after changes in 
the ionizing radiation field. 


3. NUMERICAL APPROACH 

The solution of the TDP problem is obtained by solving the 
three coupled equations (5), (13), and (16) simultaneously. To 
do so, we divide space, time, and radiation energy coordinates 
in finite elements. Thus, we express derivatives of a physical 
quantity y'd- k at the zth time step and /th spatial step as 


dy‘,j yi + lj _ yij fiyUj y< . j + 1 _ yU ) 

dt At' ’ dx Ax-i 


(24) 


with A f'7 = f' +1 7 — pT and Ax 7 = x' , -' +1 — x l, F Given the 
large temporal and spatial scales typically involved in these 
calculations, and due to the stiff nature of the differential 
equations, we find that the use of the explicit method leads to 
unstable solutions. Instead, we use the implicit method, in which 
the solution of a given equation involves both the current and a 
later state of the system. The ionization balance equation (5) is 
then expressed as: 

(ni i+1 ’') 2 [Ar i (< +1J +< +1 ’7] - n[ +1J [l +AF(2na' +1J 

+ na l c +l ’j + Y ,+l ' j )\ + \ n ‘\ + At' n 2 a l r +l '-i] = 0, (25) 
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Figure 4. As in Figure 2, but for the sudden drop in the ionizing flux by a factor of three. 
(A color version of this figure is available in the online journal.) 


where ct\ +x '^ — Q' r (7’ i+1,; ), Q!c +1 ’' , = a c (^' +1 ’■ , ), and y !+1 ' 7 = 
y(x, t l+l, j). Thus, the population n\ at the (i + l)th time 
step is given by the roots of the quadratic equation above, 
provided that the temperature r' +1 ’ 7 is known. One of these 
solutions is negative, thus non-physical, which leaves only one 
possible solution. To find the temperature we write the energy 
equation (13) as 


ji+lj 


\ + 2 a V ( ^ )2 At* 
3(n\ +1 ' J +2n‘z J ) 

2 At' 


M j 




n + 2n j 7 


3 («i +1 '' / + 2n^)k 


[n[ +1J y i+l ’j(e) 




£th] = 0. 

(26) 


The solution to this equation is found numerically by the 
secant method. Then n\ ' J is found from Equation (25) for 
every given temperature, T !+1,7 . These solutions depend on the 
photoionization rate y I+1,7 and determined through the radiative 
transfer equation (16), which in finite differences form becomes 


( C ^ t ) , T i,j, 


J 


r + J iJ ’ k (\ - cAt , n\’ J a k ) 

\2AxJ ) v 1 ' 

U+U (c&\ 

\ 2Ax> j ' 


(27) 


This equation needs to be solved for every Alh energy interval. 


Our method starts by finding the solution at t = 0, which is 
assumed to be the steady-state solution. At x = 0 the boundary 
condition is imposed: = J ‘ n k , which is the radiation field 

incident on the illuminated face of the slab. Jf£ is known at all 
times i and for every kth energy interval. 

We use logarithmically spaced grids for time, space, and 
energy. For example, for a slab of thickness Ax ~ 10 18 cm 
we use 10 3 spatial bins and a time integration over 10 4 steps up 
to t — 10 14 s, which is long enough for the system to return to 
equilibrium for all cases considered here. The resolution used 
for both the spatial and temporal grids is appropriate to resolve 
the physical phenomena relevant to this problem. We use 100 
energy bins in the 0.1-2 x 10 5 eV range. The spectral energy 
distribution of the ionizing radiation field is assumed to be a 
power law with photon index T = 2, and a high energy cut-off 
at 200 keV. 

For the present work, we investigate cases where the hydrogen 
density is kept constant at n = 10 4 cm~ 3 . Further, the intensity 
of the radiation field from the source is changed using a step 
function (i.e., instantaneous change). The change in the flux is 
specified in terms of the original flux of the source, using the 
ratio: 


fx = F™/F x , 


(28) 


where F“ ew is the new radiation flux after the change. 
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Figure 5. As in Figure 3, but for a sudden drop in the ionizing flux by a factor of three. 
(A color version of this figure is available in the online journal.) 


4. RESULTS 

4.1. Step Flux Function on a Constant Density Slab 

In this section we present simulations of photoionized slabs 

with constant hydrogen density, n = 10 4 * * * cm' 3 , subjected to a 

sudden change in the ionizing radiation. It is also assumed that 

the slabs are in steady-state equilibrium at t — 0. 

Figure 1 shows the ionization and temperature time evo- 
lutions in hydrogen clouds with two different values of §. 
The figure shows that under steady-state conditions the neutral 
hydrogen density is minimum at the illuminated side of the slab, 

where the temperature is maximum. The ionization and temper- 

ature remain relatively constant through the cloud up to a point 
where most ionizing photons have been absorbed. Then, an IF 
develops (at 7-9 x 10 16 cm) where the ionization and tempera- 
ture of the plasma drop sharply. Models for different values of f 
are very similar to each other, but the size of the ionized region 

scales up with §. 

From the state of equilibrium the incident flux suddenly 
increases by a factor of three. We follow the evolution of the 

system until it reaches a state of equilibrium again. As expected, 
a raise in the flux leads to an increment in the temperature of the 
gas and in its ionization stage, which consequently decreases 
the neutral hydrogen fraction. 

After the jump in the ionizing flux there is a temporal 
overshoot in the temperature at the IF, i.e., a sharp increase 


in temperature followed by a gradual drop to equilibrium 
values. This is due to the hardening of the ionizing flux that 
ionizes a largely neutral medium, as the lower energy ionizing 
photons get absorbed through the ionized region of the cloud. 
Moreover, a combination of fast moving photoelectrons and 
relatively few protons make recombination cooling inefficient, 
resulting in an initial sharp rise in temperature. Later, though, as 
the plasma becomes highly ionized, the recombination cooling 
rate increases, driving the temperature toward an equilibrium 
value. 

Figures 2 and 3 show the ionization/recombination rates and 
heating/cooling rates for various time steps after a change in 
the ionizing radiation field by a factor of three. At t — 0 and 
t > 3.4 x 10 8 * * * * * * * s (~10 yr) the slab is in equilibrium, thus the 
ionization and heating rates are equal to the recombination 
and cooling rates, respectively, everywhere in the cloud. In 
between these times, the figure shows ionization and heating 
fronts propagating through the cloud, leaving the plasma out of 
equilibrium. At 0.03 yrs the ionization and heating fronts are 
found at 3 x 10 16 cm, and the plasma behind these fronts is out 
of equilibrium. At t — 0.35 yrs the heating and ionization fronts 
are seen to reach the IF, where departures from equilibrium are 
at maximum. Nonetheless, by these times the gas behind the 
fronts has evolved significantly toward equilibrium. 

The ionization/recombination rates and heating/cooling rates 
are shown in Figures 4 and 5 for the case when the ionizing 
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Figure 6. Propagation time (top panel), ionization equilibration time (middle 
panel), and temperature equilibration time (bottom panel) vs. depth within the 
slab for a plasma with log £ = 0 and f x = 3. 

(A color version of this figure is available in the online journal.) 


continuum is reduced by a factor of three. Cooling and 
recombination fronts are seen to propagate through the 
cloud and behind these fronts the plasma evolves toward 
equilibrium. 


4.1.1. Timescales and Rates 

In the TDP models shown in Figure 1 the plasmas evolve 
between two steady- state solutions set by two different values 
of the ionization parameter. However, the plasma’s behavior is 
different from a sequence of equilibrium solutions calculated 
for different ionization parameters at different times. This is 
because the local conditions at different depths inside the 
cloud react at different times to the variations in the flux from 
the source, according to the propagation time. Moreover, the 
physical conditions evolve at different rates at different depths 
according to the local timescales for ionization equilibration and 
temperature equilibration. 

Figure 6 shows the propagation, ionization equilibration, and 
temperature equilibration times versus depth into the slab. It 
can be seen that fronts that result from sudden increases in 
the radiation flux travel at constant speed, ~ 20, 000 km s' 1 
(~MACH 2), from the illuminated face of the slab up to 
~3 x 10 16 cm inside the cloud. Beyond this point, the front 
slows down by orders of magnitude and the propagation time 
increases nonlinearly. In other words, it takes about ~1 yr for 
the radiation front to arrive near the IF, but several hundred years 
to move across the IF. Clearly, the absolute propagation times 
are inversely proportional to the magnitude of the flux variation, 
yet the qualitative behavior of the propagation is essentially the 
same in all cases. 

The ionization equilibration timescale depends on the relative 
change in ionization and the ionization and recombination rates. 
In steady-state conditions ionization and recombination times 
are of the order of ~100 yr, for T = 10 4 K and n e = 10 4 cm 3 . 
Thus, across the IF, where the neutral hydrogen fraction changes 
from ~1 to 0, the ionization equilibration time is about 100 yr. 
In contrast, before the IF the plasma is nearly fully ionized, thus 
the relative change in ionization is very small for any increase 
in the radiation flux and the ionization equilibration time is also 
very short. 



Distance (cm) 



Distance (cm) 


Figure 7. Ionization and temperature for a slab initially in pressure equilibrium 
at P„ = 4 x 1 0 x dyn cm' 2 . The initial flux corresponds to log $ = 0, which is 
suddenly increased by a factor of three (f x = 3). The initial condition is plotted 
in red, and the final state of the system is plotted in green. The black curves depict 
the physical conditions at different times. The gas density obtained from the 
pressure equilibrium solution is shown in the upper panel with the dashed-blue 
line. 

(A color version of this figure is available in the online journal.) 


The temperature equilibration time is of the order of a few 
years in the more ionized segment of the slab and peaks at ~35 yr 
across the IF. Interestingly, the temperature equilibration time 
is longer than the ionization equilibration time in the ionized 
fraction of the slab, but shorter across the IF. 


4.2. Step Flux Function on a Slab in Pressure Equilibrium 

Here we investigate the case of a cloud initially in gas pressure 
equilibrium with its surroundings. Let the pressure at t — 0 be 
P Q = 4 x 1(T 8 dyn cm~ 2 . For the pressure to be constant across 
the slab, the gas density increases as 1 /T from the hotter fraction 
of the cloud, facing the ionizing source, to the neutral region. 
This means that a sharp rise in density is expected across the IF, 
where the temperature drops steeply. In the present simulation 
the IF is originally found at x ~ 10 17 cm. 

In Figure 7 we show the evolution of the temperature and 
pressure when the ionizing flux is increased by a factor of three 
(f x = 3) while the gas density is kept fixed. The increase in flux 
creates an ionization and thermal front that propagates through 
the slab and heats the gas beyond the original IF. Thus, the cloud 
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Figure 8. Pressure profiles in the region where the IF is formed (black lines). The red dots indicate the position of the IF at different times. Each panel corresponds to 
a different flux variation factor f x , as indicated. 

(A color version of this figure is available in the online journal.) 


is seen to go out of pressure equilibrium, particularly across the 
original IF. As a consequence, the variation in the ionizing flux 
will induce dynamical effects in the cloud. If the thermal front 
is subsonic, the cloud will expand and the density profile of the 
gas will adjust to maintain equal pressure across the cloud and 
with its surroundings. Note that if the thermal wave is subsonic 
in the ionized region in the cloud the wave is likely to remain 
subsonic across the IF. This is because the speed of the front 
across the IF decreases roughly proportionally to T, while the 
sound speed goes as 7' l/2 . On the other hand, if the thermal 
front moves supersonically the gas has no time to adjust itself 


and strong pressure imbalances, like those seen in Figure 7, 
will appear. Thus, shocks will be formed in the slab, which can 
ultimately result in the fragmentation of the cloud (Bautista & 
Dunn 2010). Either way, variations in the ionizing flux will have 
important kinematic effects on the cloud. 

We further studied front propagations under different condi- 
tions. Figure 8 shows the pressure profiles at IFs when the flux is 
varied by factors of f x = 0.3, 0.5, 0.8, 1.2, 1.5, and 2. On each 
curve we identify the inflection point, i.e., the most negative 
value for < iP/dx, which we will use as a point of reference to 
follow up the evolution of the front. When the flux is reduced, 
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Figure 9. Propagation speed of the IFs (top panel) and recombination fronts 
(bottom panel). Each curve corresponds to a different flux variation factor/^, as 
indicated in each panel. 

(A color version of this figure is available in the online journal.) 


recombination fronts are formed and travel in the direction of 
the ionizing source. Conversely, increased fluxes lead to IFs that 
travel away from the source. 

Figure 9 shows the speeds of ionization and recombination 
fronts. It is found that the IFs move forward over long periods 
of time with speeds proportional to the flux increment (up to 
10 3 km s 1 for f x = 3). This is consistent with i> pr0 = F x /Hh 
(see Equation (23)). On the other hand, recombination fronts 
propagate with maximum speeds of the order of hundreds of 
km s' 1 for f x =0.8 or smaller. The speed of sound is given by 
v s = <Jyp/p, where y is the adiabatic index, p is the pressure, 
and p is the mass density of the gas. For an ideal gas y — 5/3 
and temperature range T = ( 1=1) x 10 4 K, v s = 12-24 km s _1 . 
Thus, even small variations of the incident flux can induce 
ionization/recombination fronts that propagate at supersonic 
speeds. 

4.3. Periodically Varying Flux on a Constant Density Slab 

In Section 4.2 we showed that equilibration times at different 
positions of a slab range by at least an order of magnitude. 
Thus, there is a large variety of astronomical nebulae whose 
radiation sources vary periodically on timescales comparable 
to their equilibration times, e.g., circumstellar nebula around 
pulsating stars and binary systems. There are also systems, like 


quasars and AGNs, characterized by quasi-periodic variability 
on all timescales. Thus, it is interesting to look at the general 
behavior of such systems. 

As discussed in previous sections, slabs with total hydrogen 
densities of ~10 4 cm - ’ have equilibration times ranging from 
less than a year to a few decades. Let us consider constant 
density slabs ionized by step-like periodically varying radiation 
continua. Figure 10 shows the neutral hydrogen density and 
temperature for various flux variation periods. This figure 
shows the average physical conditions and their full range 
of variability. For reference, we also show the steady state 
solutions for the low and high flux states and the mean conditions 
between these. Several conclusions can be drawn from this 
figure. 

1. The time average of the physical conditions is different 
from the mean of the two steady-state solutions. In general, 
the cloud tends to be overionized with respect to the steady- 
state solutions for a mean value of the flux. This is because 
ionization for a given increase in the radiation flux is a 
faster process (directly proportional to the change in the 
flux) than recombination when the flux decreases (set by 
the recombination rate coefficients and the gas density). 
On the other hand, the time-averaged temperature is lower 
than the mean of steady-state solutions in the ionized region 
of the cloud. 

2. The dispersion from the time average of the physical 
conditions increases with the period of the radiation flux. 
This is expected because for flux periods shorter than the 
plasma’s equilibration times, the cloud is forced to remain 
around a non-equilibrium state in between the two steady- 
state solutions. As the period of the flux variation increases 
the plasma has time to approach the steady-state solutions. 
Though, note that the equilibration time across the IF is 
significantly longer than the ionized region. 

3. TDP leads to much wider IFs than under steady-state 
conditions. This is due to a combination of strong gradients 
in equilibration and propagation times across the front. 
Thus, time-averaged conditions across the IF transition 
more smoothly from the ionized to neutral regions of the 
slab than under steady-state conditions. A caveat to this 
conclusion is that while the average of physical conditions 
is relatively smooth the absolute instantaneous conditions 
are not so. It is shown below that the IF exhibits larger 
variability with respect to average values than anywhere 
else in the cloud. 

Note that the behaviors discussed above are for the case of 
pure hydrogen, optically thin nebulae. Should one expect quali- 
tatively similar effects in more realistic, i.e., chemically hetero- 
geneous and optically thicker, clouds? Adding other chemical 
elements to the gas is expected to enhance cooling rates and 
optical depths. These changes are expected to have opposite ef- 
fects in terms of temporal variability. Larger cooling rates will 
contribute to reducing the temperature equilibration time. In 
turn, faster temperature equilibration will tend to drive faster 
ionization equilibration for neutral species; however, higher 
ionization stages tend to have smaller photoionization cross 
sections and for these the ionization equilibration times may 
be longer. Increasing optical depths would result in reducing 
effective recombination rates, for example by suppressing Ly« 
photons, hydrogen recombination would be reduced by ~40% 
to Case B rates, which would extend the ionization equili- 
bration times. Moreover, larger optical depths would extend 
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Figure 10. Ionization and temperature solutions for a constant density slab subjected to periodically varying fluxes with periods of 3, 9, 15, and 40 yr. The initial 
hydrogen density is 10 4 cm -3 , the radiation flux corresponds to log £ = 0, and the flux variations are of f x = ±0.5. The green curves show the steady-state equilibrium 
conditions at the low and high states of the flux. The red curves depict the steady-state equilibrium solutions for a radiation flux at the media between the low and high 
states. The blue solid line shows the time average conditions, while the dashed lines show the dispersion in that average. 

(A color version of this figure is available in the online journal.) 


propagation times in general, although the effects would vary 
along the electromagnetic spectrum and would affect different 
species selectively. In conclusion, one should expect the ef- 
fects of periodically varying continuum discussed here to be 
qualitatively valid in realistic astrophysical nebulae, albeit with 
considerable additional complexity, which deserves additional 
studies with more complete models. 

In a gas cloud photoionized by a time-dependent radiation 
source, the physical conditions change asynchronously across 
the cloud. Full animations of the ionization and temperature 
can be found at http://hea-www.cfa.harvard.edu/~javier/tdp for 
various flux variability periods. Figures 11 and 12 show a few 
snapshots of ionization and temperature conditions, normalized 
to the average values, for simulations run over 1000 yr. It is 


seen that even for radiation flux periods as long as 30 yr the 
system stays out of equilibrium through the whole duration of 
the simulation. The ionized region of the slab, which starts from 
the illuminated face, is seen to vary in sync with the continuum 
flux. On the other hand, there is a delay between the response 
across the cloud. Therefore, at any given instant, one can find, 
for example, that while most of the cloud is warmer than the time 
average, the gas across the IF would be cooler than the average. 
In general, gas across the IF behaves very differently from the 
rest of the cloud and exhibits the largest dispersion with respect 
to time-averaged conditions. This is due to the combination of 
the long propagation time and equilibration times across the IF. 
Moreover, at no time during the evolution do the gas conditions 
follow a steady-state solution. 
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Figure 11. Instantaneous ionization relative to time-averaged values for instants along 1000 yr long simulations for various radiation flux variability periods. Here, 
the radiation flux corresponds to log £ = 0 and the amplitude of variations is f x = ±0.5. 

(A color version of this figure is available in the online journal.) 
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Figure 12. Instantaneous temperature relative to time-averaged values for instants along 1000 yr long simulations for various radiation flux variability periods. Here, 
the radiation flux corresponds to log £ = 0 and the amplitude of variations is f x = ±0.5. 

(A color version of this figure is available in the online journal.) 


5. CONCLUSIONS 

We have studied the general behavior of TDP models. Here, 
the energy balance, ionization balance, and radiation transfer 
equations are considered in their full time-dependent form. 


These equations are solved for pure hydrogen plasmas subjected 
to sudden variations in the ionizing radiation field. 

Simulations of constant density slabs show the formation 
of ionization/thermal fronts that propagate through the cloud 
after a change in the ionizing flux. The propagation times and 
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response times to such fronts vary greatly from the illuminated 
face of the cloud to the IF. Simulations carried out for different 
degrees of ionization showed that the time evolution of physical 
conditions in the plasma differs from a sequence of equilibrium 
solutions. 

Our results for slabs initially in pressure equilibrium show 
that the thermal fronts that propagate through the plasma after 
a change in the ionizing flux are also pressure fronts, which 
become particularly pronounced across the IF of the slab. For 
an increase in the ionizing flux, the speed of the thermal front 
is proportional to the incident radiation flux. Thus, there is no 
limit for how fast these fronts can propagate. In contrast, a 
sudden drop in the ionizing flux creates a cooling/recombination 
front whose speed is determined by the recombination rates. 
In either case, the present simulations show that these fronts 
often propagate with supersonic speeds, and thus large pressure 
imbalances are created across the slab. This is expected to have 
important dynamical effects on the cloud, such as the creation 
of shocks and cloud fragmentation. 

Further, we studied the case of periodic variations in the 
ionizing flux. It was found that the physical conditions of the 
plasma have complex behaviors that differ from steady-state 
solutions. Moreover, even the time-averaged ionization and 
temperature are different from any steady-state case. This time 
average is characterized by overionization and a very wide IF 
with respect to the steady-state solution for a mean value of the 
radiation flux. Around the time average of physical conditions 
there is a large dispersion in the instantaneous conditions, 
particularly across the IF, which increases with the period of 
radiation flux. Moreover, the dispersion in physical conditions is 
asynchronous along the slab due to the combination of nonlinear 
propagation times for thermal/ionization front and equilibration 
times. 

Our current description of TDP is simplified owing to the 
lack of chemical elements other than hydrogen. More realistic 
models including realistic chemical mixtures and detailed mi- 
crophysics of multi-level atomic systems will be the subjects of 
future publications. 

We thank D. Hillier and P. Flarrington for useful discussions 
and valuable suggestions. 
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